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Abstract 

Using collision driven discrete molecular dynamics (DMD), we investigate the thermodynamics 
and dynamics of systems of 500 dumbbell molecules interacting by a purely repulsive ramp-like 
discretized potential, consisting of n steps of equal size. We compare the behavior of the two 
systems, with n = 18 and n = 144 steps. Each system exhibits both thermodynamic and dynamic 
anomalies, a density maximum and the translational and rotational mobilities show anomalous be- 
havior. Starting with very dense systems and decreasing the density, both mobilities first increase, 
reache a maximum, then decrease, reache a minimum, and finally increase; this behavior is similar 
to the behavior of SPC/E water. The regions in the pressure-temperature plane of translational 
and rotational mobility anomalies depend strongly on n. The product of the translational diffusion 
coefficient and the orientational correlation time increases with temperature, in contrast with the 
behavior of most liquids. 
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I. INTRODUCTION 



Recently much attention has focused on phase behavior of single component systems 
comprised of particles interacting via core-softened potentials 

Core-softened potentials exhibit a repulsive hard core plus a softening region, which can 
be a linear or nonlinear repulsive ramp, a shoulder, a double (or multiple) attractive well, 
or combinations of all these features. These models were created with the purpose of con- 
structing a simple two-body isotropic potential^ capable of describing some aspects of the 



anomalous behavior of complex fluids 
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like maximum in density clS 8b function of 



temperature, increasing isothermal com pre ssibility upon cooling, and increasing molecular 
diffusivity with the increase of pressure |ll[ . It has been proposed some time ago that these 
anomalies might be associated with a critical point at the terminus of a liquid-liquid line, 
in the unstable supercooled liquid region [3]. The study of core-softened potentials gener- 
ates models that are computationally (sometimes even analytically) tractable and that may 
retain some qualitative features of network forming fluids such as water. 

Hemmer and Stell using the method of Takahashi proposed the possibility of a 
second critical point in addition to the normal liquid-gas critical point, for a one-dimensional 
system whose pairwise potential possesses a region of negative curvature with respect to 
separation distance in its repulsive core. Debenedetti et al. developed thermodynamic 
arguments to show that this type of potential can lead to density anomaly 16]. Stillinger 
and Stillinger [l3| showed that a pure repulsive Gaussian core-softened potential produces 
both thermodynamic and dynamic anomalies but no indication of the second critical point 
at positive temperatures. 

A number of works investigated step potentials consisting of a hard core, a square re- 
pulsive shoulder and an attractive square well |2ll 122, |23|, |2|. In 
two dimensions, such potentials have density and diffusion anomalies, a negatively sloped 
freezing line and possibly a second critical point in a deeply supercooled state inaccessible by 
simulations. In three dimensions, these potentials do not have dynamic and thermodynamic 
anomalies but possess a second 13] and sometimes a third [3| critical point, accessible by 



simulations in the region predicted by the hypernetted chain integral equation |5l l25|. 

Jagla proposed a different version of the core-softened potential, the main element of 
which is a linear repulsive ramp P]. These potentials show a region of density anomaly in 
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the pressure-temperature plane and also display a liquid-liquid critical point, which can be 
located in a stable or metastable fluid region, depending on the choice of parameters 
The Rogers- Young approximation, together with the mode coupling theory, reproduces 
these features analytically [28]. An attractive ramp 3,0,131 displays a normal gas-liquid 
critical point and also brings this second critical point into an accessible region of higher 
temperature. 

An interesting issue, besides the thermodynamic anomalies, is the investigation of the 
dynamic anomalies, such as those exhibited by water. The translational diffusion coefficient 
D for a spherical particle solute of radius cr in a medium of shear viscosity r], is predicted 
by the Stokes-Einstein relation to be 



D 



(1) 



The rotational autocorrelation time is predicted by the Debye-Stokes-Einstein relation to be 

r. = ^. (2) 

Simulations of supercooled SPC/E water at constant temperature show that D increases 
as pressure P decreases, reaches a maximum at a density pomax and decreases until D reaches 



a minimum at Pn-mm 
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3^ . Further , has a minimum at pDmax and a maximum at 
and the product TrD is almost constant with density and temperature 
32| . This constancy is still an open issue and needs to be resolved. 
Combining Eq. Q with Eq. (^, we see that the product r^D is constant, which could 
explain the constancy of the product Tj.D observed the SPC/E water. However, the hydro- 
dynamic approximations used in both the Stokes-Einstein and the Debye-Stokes-Einstein 
relations are valid only if the size ratio between the solute and the solvent is large 

0- This 

is not the case of the self diffusion of SPC/E water where the tracer particle has the same 
size as the other particles in the system, so hydrodynamic arguments cannot explain the 



constancy of r^D observed in water 
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Furthermore, for a number of supercooled liquids, experimental results show that there 
is a breakdown of the Stokes-Einstein relation close to the glass transition temperature Tg 
j3J, |35|, |36|. The translational diffusion close to Tg behaves as D oc T /r]^ with 

& < 1 Q- 

n this case, the product Tr-D oc rj^~^ is not constant, but increases as the system is cooled 
3^ . Since water differs from other supercooled liquids by the presence of thermodynamic 



anomalies, the study of the relation between the thermodynamic and dynamic anomalies in 
water might shed some light on this problem. For SPC/E water, the region in the P — T 
phase diagram where water has diffusivity anomalies defines a broad region, and the region 
where density anomalies can be found lies entirely inside the anomalous diffusivity region 
. For SPC/E water, the product t^D is approximately constant inside the 
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entire region of the diffusivity anomaly. 

The goal of our study is to test the hypothesis that the constancy of TrD is the consequence 
of the diffusivity anomaly. To achieve this goal we construct a simple model which has the 
thermodynamic and dynamic anomalies of water but also has rotational degrees of freedom, 
so that Tr can be investigated. We based our model on the Jagla repulsive ramp model 
which has the entire spectrum of water anomalie s 12811 , including the correlation between the 
translational and orientation order parameters j38|. In order to create rotational degrees 
of freedom, we linked two ramp particles into a "dumbbell" by a permanent bond equal in 
length to the hard core diameter. Thus, in our model the two monomers of the dumbbell 
touch each other with their hard-cores as in rombohedral crystal obtained in 2^ near the 
high pressure boarder of the diffusivity anomaly region. We keep the same choice of the ratio 
of the soft-core to the hard-core as in j2^], since these ratio reproduces also the water-like 



correlation of the translational and orientational order maps [39|. Since hard-core models 
cannot be easily treated with a continuous molecular dynamics, we discretized the ramp 
by a step-function consisting of many small steps and use the discrete molecular dynamic 
(DMD) algorithm as in 3l- 

First, we test whether the thermodynamic and dynamic anomalies are preserved for 
dumbbells and how their region in the P — T diagram is alternated comparatively to the 
original monomeric particles. Since it it was shown that the dynamic properties depend 
significantly on the step of the discretized potential, we test if our conclusions are not an 
artifact of the discretization. Finally, in the density anomaly region, we can test if the 
product TrD is constant (as in water) or if it resembles the behavior of other supercooled 
liquids. 

This paper is organized as follows. We introduce the dumbbell model and give details 
about the simulations in Sec. II. We present the thermodynamic and dynamic behavior in 
Sec. Ill, and our conclusions in Sec. IV. 
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II. THE DUMBBELL MODEL 



The model we study is defined as follows. We consider a set of diatomic molecules 
formed by two spherical atoms of diameter (Tq linked together. The atoms in each molecule 
are separated from each other by a distance ^ which is allowed to fluctuate in the range 
0.99(To < ^ < 1.01(70. The interaction potential between two atoms belongin g to different 
dumbbell molecules separated by a distance r is a "ramp" made of n steps j28| (Fig. 1). 
Thus 



U{r) 




(3) 



where 



Ar 



(Ti - 



k - - ] Ar 



and 



AU = 



n+ (1/2) 



(4) 
(5) 

n + (l/2)- 

The liquid phase part of the ramp phase diagram has the following characteristics. At 
low P, particles prefer to be at distances Ui from each other. At high P the typical distance 
is (Jo- This implies a crossover region of P with anomalously large isothermal compressibility 
kT which might give rise to an accessible liquid-liquid phase transition if an attractive part 
is included in this potential. 



The discretized, step version of the ramp potential 
for study using collision-driven molecular dynamics 



has the advantage of being well suited 



2Ct 401 1. As n ^ oo, the step potential 

n 

becomes a "ramp" similar to the potential introduced by Hemmer and Stell [1^ and studied 
by Jagla 3l- Indeed, the properties become similar to apure ramp potential if the height 
difference between two adjacent steps is less than fc^T [2^. Moreover, we expect [s^j that 
the dynamic behavior converges to the behavior of the smooth ramp when n ^ 00 as 0(n~^), 
therefore in a slower way than the thermodynamic properties, which converge as 0(n^^). 
In order to test the influence of the discretization, we study this repulsive ramp potential 



with n = 18 and n = 144 steps and (Ji/(To = 1.74. We per 
(DMD) simulations, using a collision-driven algorithm 7, 



brm "discrete molecular dynamics" 



20 



40|, for a system comprised of 
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500 dumbbell molecules (1000 "atoms") in a cubic box, with box edge ranging from L = 13(To 
to L = 15.5cro. The dimensionless temperature T*, density p*, pressure P*, translational 
diffusion coefficient D*, and mean rotational correlational time r* are given respectively by 



and 





(7) 


* 3 

p =pao, 


(8) 


p. - P^l 


(9) 


_ D{m/U,f/^ 


(10) 


, r,(f/o/m)V2 


(11) 



where /c^ is the Boltzmann constant, m the particle mass, and Uq the high of the repulsive 
ramp. 

We calculate thermodynamic and dynamic properties using simulations of 160,000 time 
steps. We analyze the thermodynamic and dynamic behavior for dimensionless temperatures 
ranging from T* = 0.151 down to T* = 0.054 for n = 18 and from T* = 0.168 down to 
T* = 0.048 for n = 144. 



III. RESULTS 

A. Density Anomaly 

We shall first investigate if our system of diatomic molecules exhibits a density maximum 
as found for a system of spherical monoatomic particles interacting through the same po- 
tential Figure 2 plots the pressure P* against density p* for fixed values of T*. We see 
that for all the states there is no mechanical instability, and the pressure is a monotonically 
increasing function of density, 

f^Uo. (12) 



dp 



T 



implying that no phase transition is present. The isotherms cross each other, which means 
that {dp/dT)p = 0, which implies a density anomaly. 

In order to locate the points where {dP/dT)p = 0, we plot P* against T* along isochores 
(Fig. 3a, 3b and 3c), fit the results to a polynomial, and calculate the minima. The line 



-) >0. (15) 

T 



connecting these minima is the temperature of maximum density (TMD) hne, which is the 
boundary of a thermodynamically anomalous region where 

and therefore an anomalous behavior of the density (similar to water) 

(^)^>0. (M) 
Because of the Maxwell relations, Eq. ()14j) implies anomalous behavior of the entropy, 

[dp). 

We can uncover the effect of increasing n on generating the thermodynamic and dynamic 
anomalies, comparing Figs. 3a and 3b. We will see that the number of steps has almost 
no influence on the shape and location of the TMD, consistent with the fast convergence of 
pressure to the linear ramp value as n — > cxo j2^|. Therefore, the curves for the system with 
n = 144 steps are expected to be almost identical to the curves for the linear ramp j^]. 
The corresponding effect of the dumbbell shape can be estimated comparing the results of 
dumbbells with n = 18 and n = 144 with the systems of 1000 monoatomic particles with 
n = 18 (Fig. 3c) and n = 144 [2^ interacting with the same potential illustrated in Eq. (jH} 
for the same state points as the dumbbell molecule system. The systems of interacting 
dumbbells have a slightly smaller region of density (and entropy) anomaly, shifted to higher 
pressures. In these cases, for densities above 0.182 and below 0.154, no minima in the 
pressure are found. 



B. Dynamics 



We now study the mobility associated with the repulsive ramp n-step potential for n = 18 
and n = 144. We calculate the translational diffusion coefficient D using the mean-square 
displacement averaged over different initial times, 

(Ar2(t)) = ([v(to + t)-r(to)]^) (16) 

where v is the center of mass of the dumbbell. Then D is obtained from the relation 

D = \im{Ar^{t))/Qt . (17) 

t— »oo 
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Figure 4 shows the behavior of the dimensionless translational diffusion coefficient, D* ^ as 
a function of the dimensionless density, p*, for the dumbbell potential. At low temperatures, 
the behavior is similar to the behavior found for SPC/E supercooled water D increases 
as p decreases, reaches a maximum at pDmax (and Pomax); and decreases until it reaches a 
minimum at pDmin (and PDmin)- Above T* = 0.086 (not shown), however, the anomalous 
behavior disappears. 

The region in the P — T plane where there is an anomalous behavior in the diffusion 
is bounded by {Tommip), PDmm{p)) and (Tomi^Ap) , Pomi^Ap)) ■ The location of this region 
for the studied systems is shown by the dotted lines in Figs. 3a, 3b, and 3c. The effect of 
the dumbbell shape observed by comparing Fig. 3a with 3c, or Fig. 3b with Ref. [2^, is 
rather small and affects only the overall shape of the curve, shrinking slightly the region 
and shifting toward higher values of P and lower values of T. The effect of n, obtained by 

28|, however, is very stron g. W hereas in the 



comparing Fig. 3a with 3b or Fig. 3c with Ref. 
systems where the potential was discretized in 144 steps (Fig. 3b and Ref. j2S|) the diffusion 
anomaly line lies outside the TMD line (exactly as in SPC/E waterH), 

for the systems 

with n = 18 (Fig. 3a and 3c) the diffusion anomaly region shrinks strongly, in such a way 
that it migrates to a region inside the TMD line, in contrast to the behavior of SPC/E 
water. These results are consistent with the slow convergence of the dynamic properties 
of the discretized ramp potential, as n 00. Nevertheless, the results for n = 144 are 
expected to be very similar to the results for a linear ramp. 

The use of dumbbell molecules, unlike spherically symmetric models, allows us to cal- 
culate the rotational degrees of freedom and to estimate of the role of anisotropy. The 
rotational diffusion is analyzed by calculating the orientational autocorrelation function of 
the vector e{t) defining the orientation of the dumbbell molecule, 

Ce(t) = (e(t)-e(0)). (18) 

The orientational autocorrelation function depends on the density and temperature, and for 
short times we can fit Ce{t) with an exponential function |4l| : 

Ce(t) oc exp(-t/r^). (19) 

The orientational correlation times, r^, obtained from Eq. ()19j] depends on density in a 
roughly complementary way as the translational diffusion coefficient D. This behavior is 
shown in Figs. 5a (for the n = 18 system) and 5b (for the n = 144 system). 

8 



As p decreases, decreases, reaches a minimum, then increases, reac hing a maximum. 



3G 
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32j . The density 



A similar behavior was also observed for SPC/E water simulations 
region where the rotational diffusion is anomalous roughly coincides with the region of 
translational diffusion anomalies, as seen in Figs. 3a and 3b (dashed lines). Here again the 
anomaly region depends strongly on n. 

The translational and rotational diffusion for systems for the n-step repulsive ramp seem 
to behave in a complementary way, but in fact the product t^D, as shown in Fig. 6a, is not 
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. Indeed, t^D increases 



density or temperature independent as in SPC/E water 
with increasing temperature (Fig. 6b), in sharp contrast with most liquids. Inside the 
thermodynamically anomalous region the product Tj.D is also not a constant, and a careful 
examination of the very low temperature behavior (Fig. 6b) shows that the product increases 
upon cooling below T* = 0.048, in a way similar to many other systems close to the glass 
transition temperature |3Q]. For our systems we estimate the Mode Coupling Temperature 
to be 

^MCT ~ 0.044 ± 0.001 and using the Vogel-Fulcher ansatz D oc exp(yl/ (T — Tq)) we 
estimate Tq* = 0.038 ± 0.001. The glass transition temperature is between these two values. 

The explanation for the anomalous behavior of the product t^D is that the rotation 
of the dumbbell in the vicinity of other dumbbells may have a different activation energy 
than translation. This is quite plausible from simple geometrical considerations. In order 
to rotate, the dumbbell must overcome a larger potential barrier than to translate. Then 
hiTrD = {Ar — At)/kT, where both Ar and At should be functions of density. Indeed 
the Arrhenius plot of r^-D at high and moderate temperatures gives perfect straight lines 
with weak dependence on density (Fig. 6b). Figure 7 shows that Ar — At, which quantifies 
how strongly the product r^-D deviates from a constant, reaches its maximum close to the 
density corresponding to the closed packed arrangements of the rumps, below which the 
dumbbells may translate without overcoming any potential barrier. This density may be 
estimated for an arrangement of dumbbells consisting of triangular lattices of vertically 
placed dumbbells p = 1/ (crf/v^ + (j^cro\/3/2) = 0.153cro^^, where ctq is the hard core and o"i 
is considered a soft core diameter. The linearity of the Arrhenius plot breaks down close to 
the glass transition where t^D starts to increase upon cooling as in other liquids in which 
this behavior is usually associated with the breakdown of the Einstein-Stokes caused by the 



growth of dynamic heterogeneities near the glass transition 
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IV. CONCLUSIONS 



We have simulated a set of dumbbell molecules interacting through a purely repulsive 
ramp-like potential, discretized in n steps. We compare the cases of n = 18 and n = 144, 
where the n = 144 case has the behavior very similar to that if the linear ramp. We studied 
the density anomaly and anomalies in both the translational and rotational mobilities. We 
found that in both cases the density behaves anomalously in a certain range of pressures 
and temperatures. Comparison with simulations of monoatomic particles interacting with 
the same potential shows that the thermodynamic anomalies are not very sensitive to 
the number of steps used in the potential, and only weakly affected by the dumbbell shape. 
Thus the existence of a fixed link between pairs of atoms and the resulting anisotropy has 
only a small effect on the thermodynamic properties. 

The dynamic anomalies, however, are more sensitive, both to n and to the introduction 
of the dumbbell shape. The translational diffusion coefficient at constant temperature has a 
maximum and a minimum as a function of p, but the locus of anomalous behavior depends 
both on the particles' shape (monoatomic or dumbbell) as well as on the number of steps (a 
potential with 144 steps has a different range than a potential with 18 steps). The dumbbell 
shape causes the dynamically anomalous region to shrink. The decrease of n shrinks and 
shifts the dynamically anomalous region to lower temperatures. 

The region of dynamic anomalies lies approximately inside the region of thermodynamic 
anomalies for the n = 18 potential and outside for the n = 144 potential. We expect that 
systems with molecules interacting with a a smooth linear ramp potential would display a 
behavior similar to this last one. Despite of the fact that this last hierarchy of anomalous 
regions is similar to that found in simulations of water, some important differences should 
be pointed out. For dumbbell molecules, the product of the translational diffusion constant 
D and the orientational correlation time is not a constant, even at high temperatures, 
indicating that translation and rotation of dumbbells may have different activation energies. 
This is not the case in water in which the molecules are confined in symmetrical tetrahe- 
dral shells, so the translation and rotation are associated with switching the partners via 
bifurcation of the hydrogen bonds 30|. It would be interesting to compare the behavior 
of the Tj-D in the dumbbell model and other dimeric liquids, such as hydrogen peroxide. 
At low temperatures, the dumbbell model displays a sharp increase of r^-D, resembling the 
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behavior found for liquids close to Tg j36|. This behavior must be investigated by analyzing 
the "fastest" rotational and translational clusters 3, 43, 44|. 
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FIG. 1: Purely repulsive ramp-like discretized potential consisting of n steps of equal size. Step 
potential (solid line) and comparison with the linear ramp (dashed line). Here do is the diameter 
of a spheric 



1.5 



0.0 



' 1 

A 


1 1 1 

T*= 0.130 


1 1 1 1 1 1 




O 
□ 

o 


T* = 0.076 
T* = 0.059 
T* = 0.054 








1 






1 


1 , . . 





0.12 0.16 0.20 0.24 



* 

P 

FIG. 2: Dimensionless pressure, P*, as a function of dimensionless density, p*, for four different 
isotherms in the system with n = 18 steps (Fig. 1). Note that the pressure isotherms cross each 
other (dashed circle), which implies a density anomaly since dp/dT = 0. 
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FIG. 3: Dimensionless pressure against dimensionless temperature along isochores. Bold line: TMD 

line. Dotted line: boundary of the diffusivity extrema. Dashed line: boundary of the rotational 

mobility extrema. (a) Dumbbell molecules interacting via the repulsive discretized ramp p potential 

with n = 18 steps. From top to bottom, densities p* = 0.203, 0.190, 0.182, 0.175, 0.164, 0.157, 

0.148, and 0.134. In the inset it is shown the density p* = 0.164, in a scale suitable to see the 

minimum in the pressure, (b) The case of n = 144 steps. From top to botton, densities p* = 0.208, 

0.199, 0.190, 0.182, 0.175, 0.167, 0.161, 0.154, 0.148, and 0.142. (c) System of 1000 atoms (without 

the bonds), with interaction potential with n = 18 steps. From top to bottom, densities p* = 0.406, 

0,381, 0.364, 0.349, 0.328, 0.315, 0.296, 0.269, 0.254 and 0.244. 
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FIG. 4: (a) Diffusion coefficient as a function of density, along six isotlierms (from top to bottom 
T* = 0.108, 0.086, 0.075, 0.065, 0.059, 0.054) for n = 18 steps, (b) T* = 0.105, 0.084, 0.063, 0.052 
for n = 144 steps. 



17 





* — 3K — 5le- 



> > > t> — ^ — ^ — 



-* * 



0.12 



0.16 



0.20 



0.24 



FIG. 5: (a) Orientational time as a function of density for four isotherms (from bottom to top 
T* = 0.108, 0.086, 0.065, 0.054) for n = 18 steps (b) T* = 0.105, 0.084, 0.063, 0.052 for n = 144 
steps. 
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FIG. 6: (a) Product t*D* as a function of density for four isotherms (from top to bottom T* = 
0.105, 0.084, 0.063, 0.052) for n = 144 steps (b) Product t*D* as a function of the inverse of 
temperature along isochores for n = 144 steps (circles: p = 0.208, up triangles: p = 0.175, stars: p 
= 0.154, right triangles: p = 0.142) 
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FIG. 7: Difference between tfie rotational and translational activation energies, Ar and Af, calcu- 
lated as a function of density from the product TrD. 
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